Water Observations from Space (WOfS)

The DEA water observation product maps the presense of surface water from Landsat imagery. An analysis through time can map the frequency a pixel is inundated by water, which can be used to infer the temporal and spatial statistics of flood/drought events. This notebook demonstrates how to run the DEA water observation algorithm for a given area of interest

• DEA product and algorithm details
• Reference code

This notebook requires the "WOfS Environment", which needs to be installed on the first run of this notebook by executing the following cell. Once done, set the environment in the top-right corner.

# Install the WOfS environment
#  - Only installs the environment if required (run it at least once)
#  - Select the "WOfS Environment" kernel as per the picture above
!../tools/install_wofs.sh

WOfS Environment kernel already installed

Table of contents

• Data parameters - Lake Tempe
• Set proxy parameters to access data locally
• Import dependencies and initialise datacube
• Display the region of interest
• Load and display the DEM
• Load the data
• Classify WOfS
• Water Observations Summaries
  • Wet counts
  • Clear observations counts
  • Wet frequencies

Data parameters - Lake Tempe

aoi_label = "Lake Tempe, Indonesia"
product = "landsat8_c2l2_sr"
longitude = (119.8242517, 120.0350519)
latitude = (-4.2013799, -3.9445384)
time = ('2020-01-01', '2020-12-31')
output_crs = "EPSG:32650"
resolution = (30, -30)

# Where to save the DEM fetched in ODC
DEM_PATH = "/home/jovyan/dems/srtm_lake_tempe.tif"

Set proxy parameters to access data locally

from os import environ

environ["AWS_HTTPS"] = "NO"
environ["GDAL_HTTP_PROXY"] = "easi-caching-proxy.caching-proxy:80"
print(f'Will use caching proxy at: {environ.get("GDAL_HTTP_PROXY")}')

Will use caching proxy at: easi-caching-proxy.caching-proxy:80

Import dependencies and initialise datacube

import sys
from pathlib import Path

import xarray as xr
import rioxarray

from wofs.virtualproduct import WOfSClassifier

try:
    from dea_tools.plotting import display_map, rgb
except ModuleNotFoundError:
    # Local copy of selected dea_tools
    if '../tools' not in sys.path:
        sys.path.append('../tools')
    from datacube_utils import display_map
    rgb = None  # Not copied or adapted yet

# Ignore warnings in output
import warnings
from sqlalchemy.exc import SAWarning
warnings.filterwarnings("ignore", category=FutureWarning) 
warnings.filterwarnings("ignore", category=SAWarning) 

from os import environ
from cartopy.crs import PlateCarree
from datashader import reductions
import holoviews as hv
import hvplot.xarray
import matplotlib.pyplot as plt

from datacube import Datacube
from datacube.utils.aws import configure_s3_access

configure_s3_access(
    aws_unsigned=False, 
    requester_pays=True, 
);

dc = Datacube()

# Display available products
products_info = dc.list_products()
products_info

	name	description	license	default_crs	default_resolution
name
copernicus_dem_30	copernicus_dem_30	Copernicus 30m Digital Elevation Model (GLO-30)	None	None	None
landsat5_c2l2_sr	landsat5_c2l2_sr	Landsat 5 Collection 2 Level-2 Surface Reflect...	None	None	None
landsat5_c2l2_st	landsat5_c2l2_st	Landsat 5 Collection 2 Level-2 UTM Surface Tem...	CC-BY-4.0	None	None
landsat7_c2l2_sr	landsat7_c2l2_sr	Landsat 7 USGS Collection 2 Surface Reflectanc...	None	None	None
landsat7_c2l2_st	landsat7_c2l2_st	Landsat 7 Collection 2 Level-2 UTM Surface Tem...	CC-BY-4.0	None	None
landsat8_c2l1	landsat8_c2l1	Landsat 8 Collection 2 Level-1 (top of atmosph...	CC-BY-4.0	None	None
landsat8_c2l2_sr	landsat8_c2l2_sr	Landsat 8 Collection 2 Surface Reflectance, pr...	None	None	None
landsat8_c2l2_st	landsat8_c2l2_st	Landsat 8 Collection 2 Level-2 UTM Surface Tem...	CC-BY-4.0	None	None
landsat9_c2l1	landsat9_c2l1	Landsat 9 Collection 2 Level-1 (top of atmosph...	CC-BY-4.0	None	None
landsat9_c2l2_sr	landsat9_c2l2_sr	Landsat 9 Collection 2 Surface Reflectance, pr...	CC-BY-4.0	None	None
landsat9_c2l2_st	landsat9_c2l2_st	Landsat 9 Collection 2 Level-2 UTM Surface Tem...	CC-BY-4.0	None	None
lpdaac_nasadem	lpdaac_nasadem	NASADEM Merged DEM Global 1 arc second V001	None	None	None
nasa_aqua_l2_oc	nasa_aqua_l2_oc	NASA MODIS-Aqua L2 Ocean Color, regridded to W...	None	EPSG:4326	(0.01, 0.01)
nasa_aqua_l2_sst	nasa_aqua_l2_sst	NASA MODIS-Aqua L2 Sea Surface Temperature, re...	None	EPSG:4326	(0.01, 0.01)
s2_l2a	s2_l2a	Sentinel-2a and Sentinel-2b imagery, processed...	None	None	None

Display the region of interest

display_map(x=longitude, y=latitude)

Make this Notebook Trusted to load map: File -> Trust Notebook

Load and display the DEM

Load the SRTM DEM data from ODC. If the area of interest is reasonably flat (like Lake Tempe) the DEM can be optional.

In this notebook, the DEM data is retrieved and saved to the same local file each time the next cells are run. This if fine if the area of interest is being changed (as per this demonstration notebook). For other use-cases, consider retrieving and saving the DEM data once only.

dem = dc.load(
    product="lpdaac_nasadem", 
    latitude=latitude,
    longitude=longitude,
    output_crs="epsg:4326", 
    resolution=(-1/3600, 1/3600),
)
elevation = dem.elevation

options = {
    'title': 'Elevation',
    'width': 800,
    'height': 500,
    'aspect': 'equal',
    'cmap': plt.cm.terrain,
    'clim': (0, elevation.max().values.item()),    # Limit the color range depending on the layer_name
    'colorbar': True,
    'tools': ['hover'],
}
plot_crs = PlateCarree()
elevation.hvplot.image(
     x = 'longitude', y = 'latitude',         # Dataset x,y dimension names 
     crs = plot_crs,
     rasterize = True,                        # If False, data will not be reduced. This is slow to load but all data is loaded.
     aggregator = reductions.mean(),          # Datashader calculates the mean value for reductions (also first, min, max, las, std, mode)
     precompute = True,                       # Datashader precomputes what it can
    ).opts(**options).hist(bin_range = options['clim'])

dem_path = Path(DEM_PATH)
dem_path.parent.mkdir(parents=True, exist_ok=True)
elevation.rio.to_raster(dem_path)

Load the data

Load the data from ODC and rename bands as needed by the WOfS classifier.

# Ignore SAWarning in output
measurements = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
data = dc.load(
    product=product,
    longitude=longitude,
    latitude=latitude,
    time=time,
    output_crs=output_crs,
    resolution=resolution,
    measurements=measurements,
    dask_chunks={'time': 1},
)

data = data.rename({
    "blue": "nbart_blue",
    "green": "nbart_green",
    "red": "nbart_red",
    "nir": "nbart_nir",
    "swir1": "nbart_swir_1",
    "swir2": "nbart_swir_2",
    "pixel_qa": "fmask",
})
data

<xarray.Dataset>
Dimensions:       (time: 19, y: 951, x: 785)
Coordinates:
  * time          (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-1...
  * y             (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
  * x             (x) float64 8.371e+05 8.37e+05 ... 8.136e+05 8.136e+05
    spatial_ref   int32 32650
Data variables:
    nbart_blue    (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_green   (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_red     (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_nir     (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_swir_1  (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_swir_2  (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    fmask         (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
    crs:           EPSG:32650
    grid_mapping:  spatial_ref

xarray.Dataset

• Dimensions:
  • time: 19
  • y: 951
  • x: 785
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2020-01-09T02:10:26.574371 ... 2...

    units :

    seconds since 1970-01-01 00:00:00

    array(['2020-01-09T02:10:26.574371000', '2020-01-25T02:10:22.588145000',
           '2020-02-10T02:10:17.788921000', '2020-02-26T02:10:13.879802000',
           '2020-03-13T02:10:07.386121000', '2020-03-29T02:09:58.450496000',
           '2020-04-14T02:09:51.264053000', '2020-04-30T02:09:43.329680000',
           '2020-05-16T02:09:40.092585000', '2020-06-01T02:09:44.725513000',
           '2020-06-17T02:09:54.647535000', '2020-07-03T02:10:02.482792000',
           '2020-08-20T02:10:18.209098000', '2020-09-05T02:10:25.859941000',
           '2020-09-21T02:10:31.337497000', '2020-10-23T02:10:34.865932000',
           '2020-11-24T02:10:34.546460000', '2020-12-10T02:10:36.006848000',
           '2020-12-26T02:10:33.500205000'], dtype='datetime64[ns]')

  • y
    (y)
    float64
    -4.65e+05 -4.65e+05 ... -4.365e+05

    units :

    metre

    resolution :

    30.0

    crs :

    EPSG:32650

    array([-465045., -465015., -464985., ..., -436605., -436575., -436545.])

  • x
    (x)
    float64
    8.371e+05 8.37e+05 ... 8.136e+05

    units :

    metre

    resolution :

    -30.0

    crs :

    EPSG:32650

    array([837075., 837045., 837015., ..., 813615., 813585., 813555.])

  • spatial_ref
    ()
    int32
    32650

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]

    grid_mapping_name :

    transverse_mercator

    array(32650, dtype=int32)

• Data variables: (7)
  • nbart_blue
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    units :

    reflectance

    nodata :

    0

    crs :

    EPSG:32650

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 27.05 MiB 1.42 MiB Shape (19, 951, 785) (1, 951, 785) Count 38 Tasks 19 Chunks Type uint16 numpy.ndarray

  • nbart_green
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    units :

    reflectance

    nodata :

    0

    crs :

    EPSG:32650

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 27.05 MiB 1.42 MiB Shape (19, 951, 785) (1, 951, 785) Count 38 Tasks 19 Chunks Type uint16 numpy.ndarray

  • nbart_red
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    units :

    reflectance

    nodata :

    0

    crs :

    EPSG:32650

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 27.05 MiB 1.42 MiB Shape (19, 951, 785) (1, 951, 785) Count 38 Tasks 19 Chunks Type uint16 numpy.ndarray

  • nbart_nir
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    units :

    reflectance

    nodata :

    0

    crs :

    EPSG:32650

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 27.05 MiB 1.42 MiB Shape (19, 951, 785) (1, 951, 785) Count 38 Tasks 19 Chunks Type uint16 numpy.ndarray

  • nbart_swir_1
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    units :

    reflectance

    nodata :

    0

    crs :

    EPSG:32650

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 27.05 MiB 1.42 MiB Shape (19, 951, 785) (1, 951, 785) Count 38 Tasks 19 Chunks Type uint16 numpy.ndarray

  • nbart_swir_2
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    units :

    reflectance

    nodata :

    0

    crs :

    EPSG:32650

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 27.05 MiB 1.42 MiB Shape (19, 951, 785) (1, 951, 785) Count 38 Tasks 19 Chunks Type uint16 numpy.ndarray

  • fmask
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    units :

    bit_index

    nodata :

    1

    flags_definition :

    {'snow': {'bits': 5, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'clear': {'bits': 6, 'values': {'0': 'not_clear', '1': 'clear'}}, 'cloud': {'bits': 3, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}, 'cirrus': {'bits': 2, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_pixel': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 'values': {'1': 'Fill', '2': 'Dilated Cloud', '4': 'Cirrus', '8': 'Cloud', '16': 'Cloud Shadow', '32': 'Snow', '64': 'Clear', '128': 'Water', '256': 'Cloud Confidence low bit', '512': 'Cloud Confidence high bit', '1024': 'Cloud Shadow Confidence low bit', '2048': 'Cloud Shadow Confidence high bit', '4096': 'Snow Ice Confidence low bit', '8192': 'Snow Ice Confidence high bit', '16384': 'Cirrus Confidence low bit', '32768': 'Cirrus Confidence high bit'}, 'description': 'Level 2 pixel quality'}, 'cloud_shadow': {'bits': 4, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}}, 'cloud_confidence': {'bits': [8, 9], 'values': {'0': 'none', '1': 'low', '2': 'medium', '3': 'high'}}, 'cirrus_confidence': {'bits': [14, 15], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'snow_ice_confidence': {'bits': [12, 13], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'cloud_shadow_confidence': {'bits': [10, 11], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}}

    crs :

    EPSG:32650

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 27.05 MiB 1.42 MiB Shape (19, 951, 785) (1, 951, 785) Count 38 Tasks 19 Chunks Type uint16 numpy.ndarray

• Attributes: (2)

  crs :

  EPSG:32650

  grid_mapping :

  spatial_ref

Classify WOfS

Run the water observations classifier. This may take a few minutes.

transform = WOfSClassifier(c2_scaling=True, dsm_path=DEM_PATH)
# Compute the WOFS layer
wofl = transform.compute(data)
wofl

<xarray.Dataset>
Dimensions:      (y: 951, x: 785, time: 19)
Coordinates:
  * y            (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
  * x            (x) float64 8.371e+05 8.37e+05 8.37e+05 ... 8.136e+05 8.136e+05
  * time         (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-12...
    spatial_ref  int32 32650
Data variables:
    water        (time, y, x) uint8 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
    crs:      EPSG:32650

xarray.Dataset

• Dimensions:
  • y: 951
  • x: 785
  • time: 19
• Coordinates: (4)
  • y
    (y)
    float64
    -4.65e+05 -4.65e+05 ... -4.365e+05

    units :

    metre

    resolution :

    30.0

    crs :

    EPSG:32650

    array([-465045., -465015., -464985., ..., -436605., -436575., -436545.])

  • x
    (x)
    float64
    8.371e+05 8.37e+05 ... 8.136e+05

    units :

    metre

    resolution :

    -30.0

    crs :

    EPSG:32650

    array([837075., 837045., 837015., ..., 813615., 813585., 813555.])

  • time
    (time)
    datetime64[ns]
    2020-01-09T02:10:26.574371 ... 2...

    units :

    seconds since 1970-01-01 00:00:00

    array(['2020-01-09T02:10:26.574371000', '2020-01-25T02:10:22.588145000',
           '2020-02-10T02:10:17.788921000', '2020-02-26T02:10:13.879802000',
           '2020-03-13T02:10:07.386121000', '2020-03-29T02:09:58.450496000',
           '2020-04-14T02:09:51.264053000', '2020-04-30T02:09:43.329680000',
           '2020-05-16T02:09:40.092585000', '2020-06-01T02:09:44.725513000',
           '2020-06-17T02:09:54.647535000', '2020-07-03T02:10:02.482792000',
           '2020-08-20T02:10:18.209098000', '2020-09-05T02:10:25.859941000',
           '2020-09-21T02:10:31.337497000', '2020-10-23T02:10:34.865932000',
           '2020-11-24T02:10:34.546460000', '2020-12-10T02:10:36.006848000',
           '2020-12-26T02:10:33.500205000'], dtype='datetime64[ns]')

  • spatial_ref
    ()
    int32
    32650

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]

    grid_mapping_name :

    transverse_mercator

    array(32650, dtype=int32)

• Data variables: (1)
  • water
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 951, 785), meta=np.ndarray>

    Array Chunk Bytes 13.53 MiB 729.04 kiB Shape (19, 951, 785) (1, 951, 785) Count 2469 Tasks 19 Chunks Type uint8 numpy.ndarray

• Attributes: (1)

  crs :

  EPSG:32650

# Uncomment the following line to display WOfS for each time slice. This may take a few minutes
# wofl.water.plot(col="time", col_wrap=5);

Water observations summaries

Water observations summaries, based on odc-stats/plugins/wofs.py, consist of:

• count_clear: a count of every time a pixel was observed (not obscured by terrain or clouds)
• count_wet: a count of every time a pixel was observed and wet
• frequency: what fraction of time (wet/clear) was the pixel wet

# Rename dimensions as required
wofl = wofl.rename({"x": "longitude", "y": "latitude"})

from odc.algo import safe_div, apply_numexpr, keep_good_only

wofl["bad"] = (wofl.water & 0b0111_1110) > 0
wofl["some"] = apply_numexpr("((water<<30)>>30)==0", wofl, name="some")
wofl["dry"] = wofl.water == 0
wofl["wet"] = wofl.water == 128
wofl = wofl.drop_vars("water")
for dv in wofl.data_vars.values():
    dv.attrs.pop("nodata", None)

# Ignore warnings triggered by time slices without data at all
warnings.filterwarnings("ignore", message="divide by zero encountered in true_divide") 
warnings.filterwarnings("ignore", message="invalid value encountered in true_divide") 

wofl.wet.plot(col="time", col_wrap=5);

# Helper frunction from https://github.com/opendatacube/odc-stats/blob/develop/odc/stats/plugins/wofs.py
def reduce(xx: xr.Dataset) -> xr.Dataset:
    nodata = -999
    count_some = xx.some.sum(axis=0, dtype="int16")
    count_wet = xx.wet.sum(axis=0, dtype="int16")
    count_dry = xx.dry.sum(axis=0, dtype="int16")
    count_clear = count_wet + count_dry
    frequency = safe_div(count_wet, count_clear, dtype="float32")

    count_wet.attrs["nodata"] = nodata
    count_clear.attrs["nodata"] = nodata

    is_ok = count_some > 0
    count_wet = keep_good_only(count_wet, is_ok)
    count_clear = keep_good_only(count_clear, is_ok)

    return xr.Dataset(
        dict(
            count_wet=count_wet,
            count_clear=count_clear,
            frequency=frequency,
        )
    )

summary = reduce(wofl)

Wet counts

A count of every time a pixel was observed and wet.

summary.count_wet.plot(size=10)
plt.title(f"{aoi_label} – Wet observations counts for {time[0]} to {time[1]}");

Clear observations counts

A count of every time a pixel was observed (not obscured by terrain or clouds).

summary.count_clear.plot(size=10)
plt.title(f"{aoi_label} – Clear observations counts for {time[0]} to {time[1]}");

Wet frequencies

What fraction of the time was the pixel wet.

summary.frequency.plot(size=10)
plt.title(f"{aoi_label} - Wet frequency for {time[0]} to {time[1]}");

# Save frequencies to high-res image file
summary.frequency.plot(size=20)
plt.title(f"{aoi_label} - Wet frequency for {time[0]} to {time[1]}")
plt.savefig('/home/jovyan/tempe_wofs_2020.png', dpi=600);